clear all
clc

load('D:\My Documents\CIRANO\db_mySQL\tp_2011Dec11_MR7.mat','conditions','numEdges','year_list')

conn = database('BDC','','');
%yr = 11; %2008
yr = 9;%2006
%yr = 8;%2005
%yr = 4;%2001
%yr = 5;%2002
%yr = 3;% 2000

T = 1; % time horizon
R = 0.45; % recovery rate
numSimulate = 5000;
SG = cellstr(char('Sample 1: [$50k,+)','Sample 2A: [$50k, 250k)','Sample 2B: [$250k, 1000k)','Sample 2C: [$1000k, +)'));

loss = zeros(numSimulate,1);
transMatGrp = zeros(length(conditions),numEdges,numEdges);
portfolioGrp = zeros(length(conditions), numEdges-1, 3);
defaultGrp = zeros(length(conditions), numSimulate,numEdges-1);
lossGrp = zeros(numSimulate,length(conditions)); 
transMatGrp(1,:,:) = xlsread('D:\My Documents\R\RResult\Q1_ml_06_WRcensor_4grps_8states.csv','A2:H9');
transMatGrp(2,:,:) = xlsread('D:\My Documents\R\RResult\Q2_ml_06_WRcensor_4grps_8states.csv','A2:H9');
transMatGrp(3,:,:) = xlsread('D:\My Documents\R\RResult\Q3_ml_06_WRcensor_4grps_8states.csv','A2:H9');
transMatGrp(4,:,:) = xlsread('D:\My Documents\R\RResult\Q4_ml_06_WRcensor_4grps_8states.csv','A2:H9');
%transMatGrp(4,:,:) = xlsread('D:\My Documents\R\RResult\Q4_ml_08_deathCensor_4grps.csv','A2:H9');
%transMatGrp(4,:,:) = xlsread('D:\My Documents\R\RResult\Q4_ml_08_4grps.csv','A2:H9');
for grp = 1:length(conditions)
    
% calcuate exposure of assigned year
condition = char(conditions(grp));
sqlquery1 = strcat(['select grpRate(modiRate) as grpRate, avg(PPDOS_sum), count(customer) from helena20111207.final_v4_tbl as t where Year(modiFisDate) =  ''' num2str(year_list(yr)) ''' and Month(modiFisDate) = 12 and ' condition ' group by grpRate;']);
exposure = fetch(conn, sqlquery1);

% set portfolio distribution of rating and exposure
portfolio (:,:) = cell2mat(exposure(2:numEdges, :));
portfolioGrp(grp, :,:) = portfolio(:,:);
display(['Portfolio total firm# in group ' num2str(grp) ' is ' num2str(sum(portfolio(:,3)))])

% get generator for grp and yr
diagInd = logical(eye(numEdges,numEdges));
transMat(:,:) = transMatGrp(grp,:,:);

% % set subremoveDs--This maybe wrong, should not sort--Mar1,2012
% subremoveD = bsxfun(@rdivide, transMat, -transMat(diagInd));
% [vs, vi] = sort(subremoveD,2); % vs is sorted vectors, vi is indexed vectors

% set subremoveDs_Mar/1/2012
subinterval = bsxfun(@rdivide, transMat, -transMat(diagInd));
removeD = removeDigonal(subinterval);
vs = zeros(size(removeD));
vs(:,1) = removeD(:,1); % vs is subinterval vector
for i = 2: size(removeD,2)
vs (:,i) = vs(:,i-1)+removeD(:,i);
end
vi = repmat([1:8],8,1); % vi is index of vs
vi = removeDigonal(vi);


% simulation
display(['start simulation for group ' num2str(grp)]);

for i = 1: numSimulate
    if mod(i,100) == 0
        disp( strcat(['Group ' num2str(grp) ': i= ' num2str(i)]));
    end
    default = zeros(numEdges-1,1);
    for r = 1: numEdges-1
        for f = 1: portfolio(r,3)
            t = 0;
            startR = r;
            t = t+ exprnd(-1/transMat(startR,startR));
            
            while t < T
                x = unifrnd(0,1);
                jumpR = jumpRate(x, vs, vi, startR);
                if ~(jumpR == 8)
                    startR = jumpR;
                    t = t+ exprnd(-1/transMat(startR,startR));
                else
                    default(r,1) = default(r,1)+1;
                    %Loss_temp = [Loss, portfolio(r,2)* (1-R)];
                    break; %Mar/1/2012 added to correct the default case 
                end
            end
        end
    end

defaultGrp(grp, i, :) = default;
% loss(i,1) = sum(default.* portfolio(:,2)*(1-R)); % correct it already! before i wrongly put portfolio(:,3) which is not exposure
% lossGrp(i,grp) = loss(i,:);
end

% 
% subplot(2,2,grp);hist(loss(:,1));
% title(SG(grp));
% xlabel('Loss');
% ylabel('Frequency')    
% grid on

end

close(conn);

%% To correct the wrong result for loss distribution in previous coding

loss_right = zeros(numSimulate,1);% for correction in 2nd section
lossGrp_right = zeros(numSimulate,length(conditions)); % for correction in 2nd section
for grp = 1: 4
    for i = 1: numSimulate
        default_temp(:,:) = defaultGrp(grp, i, :);
        portfolio(:,:)= portfolioGrp(grp, :,:);
        loss_right(i,1) = sum(default_temp.* portfolio(:,2)*(1-R));
        lossGrp_right(i,grp) = loss_right(i,:);
    end
end

miu = mean(lossGrp_right);
sig =  std(lossGrp_right,0,1);
VaR = quantile (lossGrp_right,[.99 .95]);
%VaR = prctile(lossGrp_right,[1 5]); <--wrong! should be quantile and 99&95

% for grp = 4:4
% subplot(2,2,grp);
% a = lossGrp_right(:,grp);
% hist(a);
% title(SG(grp));
% xlabel('Loss');
% ylabel('Frequency')    
% grid on
% end
save('D:\My Documents\CIRANO\db_mySQL\simulation_06_5k_ml_correctMar2012_WRCensor_8states_yes.mat')
%save ('D:\My Documents\CIRANO\db_mySQL\simulation_08_5k_ml_correct_grp4_deathCensor.mat')
%save ('D:\My Documents\CIRANO\db_mySQL\simulation_08_5k_ml_correct_grp4.mat')



% [y,xi] = ksdensity(loss);
% plot(xi,y); 